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Radiative  Transfer  in  Falling  Snow 

A  Two-stream  Approximation 


GARY  KOH 


INTRODUCTION 

The  effects  of  multiple  scattering  must  be  considered  to  describe  the  transmission  of  visi¬ 
ble  light  through  an  optically  thick  atmosphere  containing  haze,  fog,  clouds  or  solid  aerosol 
particles.  The  rays  or  the  photons  that  propagate  through  the  atmosphere  may  undergo 
many  scattering  events  before  entering  or  escaping  the  field  of  view  of  a  detector  located 
downstream.  Simple  analytical  expressions  that  are  easy  to  interpret  and  that  adequately  rep¬ 
resent  the  important  features  of  the  multiple  scattering  process  can  be  obtained  using  various 
two-stream  approximations  to  radiative  transfer.  Typically,  the  two-stream  models  param¬ 
eterize  the  radiative  properties  of  the  atmosphere  so  that  the  intensities  of  the  reflected  and 
transmitted  radiation  can  be  calculated  as  a  function  of  the  optical  depth,  the  single  scatter¬ 
ing  albedo  and  the  moments  of  the  phase  function. 

The  two-stream  approximations  are  usually  derived  by  beginning  with  Chandrasekhar’s 
(1960)  exact  equation  of  radiative  transfer  and  by  making  various  approximations  (Meador 
and  Weaver  1980,  King  and  Harshvardhan  1986).  Bohren  (1987)  presents  a  derivation  of  a 
simple  two-stream  approximation  in  which  the  photons  are  constrained  to  be  scattered  strict¬ 
ly  in  the  forward  and  backward  directions.  He  avoids  the  use  of  Chandrasekhar’s  equation 
and  gives  a  derivation  similar  to  that  of  Schuster  (1905),  who  initiated  the  two-stream  radia¬ 
tion  transfer  theory  long  before  the  exact  equation  was  developed.  This  report  derives  the 
two-stream  approximation  presented  by  Bohren  (1987)  by  starting  with  the  exact  equation  of 
radiative  transfer  to  show  the  relationship  between  the  approximate  and  the  exact  solutions. 
The  simple  two-stream  theory  is  then  used  to  describe  the  effects  of  multiple  scattering  on  a 
parallel  beam  of  light  transmitted  through  falling  snow.  The  comparison  of  the  approximate 
results  with  the  experimental  results  shows  that  the  simple  two-stream  theory  can  be  used  to 
describe  the  transmission  of  visible  light  through  falling  snow. 


DESCRIPTION  OF  THE  PROBLEM 

A  parallel  beam  of  light  incident  on  a  snow-filled  atmosphere  is  scattered  and  absorbed  so 
that  the  intensity  of  light  reaching  the  detector  located  downstream  decreases.  To  describe 
this  decrease  in  intensity,  it  is  convenient  to  distinguish  between  the  incident  intensity,  /j,  the 
reduced  (unscattered  and  unabsorbed)  intensity,  4,  the  diffuse  (scattered)  intensity,  /j,  and 
the  total  intensity,  /(.  Figure  1  is  a  diagram  illustrating  the  problem. 


X=0  'C  =  0«rt-L 


Figure  I.  Relationship  between  the  in¬ 
cident  intensity  (Ij),  the  reduced  inten¬ 
sity  (If)  and  the  diffuse  intensity  (i^.  r 
is  the  optical  depth,  L  is  the  path  length 
and  fi  =  cos$. 


TliC  reduced  intensity  is  the  portion  of  the  incident  light  that  travels  the  distance  L  without 
undergoing  scattering  and  absorption  and  satisfies  the  equation 

dl,  =  (1) 

where  is  the  extinction  coefficient.  The  extinction  coefficient  is  defined  as  the  sum  of  the 
absorption  coefficient  and  the  scattering  coefficient.  This  equation  integrates  to  yield  the  fa¬ 
miliar  Beer-Lambert  law  and  the  reduced  intensity  is  expressed  as 


I,  =  /j  exp(  -  t) 


(2) 


where  the  optical  depth,  t,  is  defined  as 

The  diffuse  intensity  represents  the  portion  of  the  incident  light  that  is  scattered  one  or 
more  times  and  that  travels  in  all  directions  according  to  the  scattering  characteristics  of  the 
snow  particles.  The  equation  that  describes  the  transmission  of  diffuse  intensity  through  a 
plane  parallel  atmosphere  is  the  integro-differential  equation  of  radiative  transfer  (Chand¬ 
rasekhar  1960,  p.  22) 

d  1  * 

-/j  exp(-T/^io)  /i,0)  (3) 

where  /<  is  the  cosine  of  the  angle  0  measured  from  the  axis  normal  to  the  parallel  atmos¬ 
phere,  and  o  is  the  corresponding  azimuthal  angle.  The  phase  function  e(/i',0':  de¬ 

scribes  a  single  scattering  event  from  the  direction  ',0 '  into  the  direction  n,()>.  The  incident 
beam  comes  from  -ft,  direction,. 

The  total  intensity  is  the  sum  of  the  reduced  intensity  and  the  diffuse  intensity  and  is  ex¬ 
pressed  as 


A  =  A  +  ^d-  (4) 

From  eq  2  and  3  it  is  seen  that  the  reduced  intensity  diminishes  exponentially  while  the  be¬ 
havior  of  the  diffuse  intensity  is  more  complex.  To  avoid  the  extensive  computations  that  are 
required  to  solve  for  the  diffuse  intensity,  a  simple  two-stream  approximation  is  used. 


TWO-STREAM  APPROXIMATION 

The  two-stream  approximation  used  in  this  study  is  obtained  by  making  several  assump¬ 
tions.  It  is  initially  assumed  that  the  diffuse  intensity  is  azimuthally  independent  and  con¬ 
stant  over  the  forward  and  the  backward  hemisphere  (Schuster’s  approximation)  so  that 

h  >  0 

[d{T,ti,6)  =  (5) 

<  0 

where  1^(7, ti)  and  ld(T,n)  represent  the  intensities  for  the  diffuse  radiation  in  the  forward 
and  backward  hemispheres  respectively.  The  diffuse  radiation  is  then  constrained  to  travel  in 
two  directions  only,  exactly  forward  and  exactly  backward,  so  that  the  intensities  can  be  ex¬ 
pressed  as  /(f(T,  -t- 1)  and  ld(T,  - 1)  respectively.  Applying  these  constraints  to  eq  3  results 
in  the  following  two  coupled  linear  differential  equations 


2 


^ /d  (’■)  =  - +^( “  +  ^) ^d (’■)+^(  + 1;  + 1) j ~ 


+  F(-l;  +  l)/iexp(-T) 


(6a) 


£  liir)  =  /d'^Cr)  -  F(  + 1 :  -  1)  /d"(T)  -  F(  - 1 :  -  1)  /j-'Cr) ;  m  =  "  1 


-F(-l:-l)/iexp(-T)  (6b) 

where  I^t)  represents  the  light  that  comes  from  -  1  direction  and  travels  toward  +  1  direc¬ 
tion  and  I^t)  represents  the  light  that  comes  from  -t- 1  direction  and  travels  toward  -  1  di¬ 
rection.  Fifi  '\ii)  is  the  fraction  of  light  incident  on  a  snow  particle  from  the  n  '  direction  that 
is  scattered  in  the  /t  direction.  For  example,  F(-  1,  -  1)  and  F(-(- 1,  -(-1)  represent  the  frac¬ 
tion  of  the  light  that  is  scattered  in  the  direction  from  which  the  light  came  (backscattering). 

The  values  of  Fiji  ';^)  approximate  the  integration  in  eq  3  for  the  forward  and  backward 
hemispheres  so  that 


F(-l;  +  l) 

F(+l;  +  l) 


A  1  f  Qin',4>'\ii,4>)d<t>'dn' 


4ir 


-1  0 


$  \  Qil^' ,<t>'',lJ;<t>)d<t>'dn' 


Air 


0  0 


(7a) 

(7b) 


F(-t- 1;  -  1)  and  F(+  l;  +  I )  can  be  similarly  represented. 

The  phase  function  of  snow  is  not  known,  therefore  the  integration  in  eq  7a  and  7b  must 
be  carried  out  using  some  approximate  phase  function.  The  merits  of  performing  such  com¬ 
plex  integrations  using  approximate  values  for  the  phase  function  are  questionable.  F(/t';/i) 
may  be  more  conveniently  approximated  in  terms  of  the  asymmetry  parameter,  g,  of  the 
phase  function  and  the  single  scattering  albedo,  wo,  so  that  (van  de  Hulst  1980) 


F(  +  l;-l)  =  F(-l;  +  l)  =  (8a) 

F(-t-l;  +  l)  =  F(-l;-l)  =  (8b) 

The  asymmetry  parameter,  g,  is  the  mean  cosine  of  the  scattering  angle  and  defines  the 
degree  of  anisotropy  of  a  phase  function.  The  values  for  g  range  from  -  1  to  -I- 1 . 

For  isotropic  scattering  and  symmetrical  scattering  about  90°,  the  asymmetry  parameter  is 
0  so  that  the  values  of  F(^  ';/i)  in  eq  8a  and  8b  are  one-half.  In  other  words,  an  isotropic  scat- 
terer  scatters  one-half  of  the  incident  light  in  one  direction  and  the  other  half  in  the  opposite 
direction.  For  a  highly  anisotropic  scatterer  in  the  forward  direction  the  value  for  g  ap¬ 
proaches  1 .  In  this  case  eq  8a  and  8b  approach  1  and  0  respectively.  This  means  that  most  of 
the  light  is  scattered  in  the  direction  of  the  incident  beam  (forward  scattering).  Similar  rea¬ 
soning  can  be  applied  to  a  scatterer  whose  asymmetry  parameter  approaches  -  1 . 

By  substitution  of  the  results  from  eq  8a  and  8b  into  eq  6a  and  6b,  a  two-stream  model  can 
be  expressed  by 


£  /d(r)  =  -  /d(r)  +  w.  /d-(r)  +  w,  /j  exp(  -  t)  (9a) 
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~  /<rW  =  /d-W-oJo  (9b) 

These  coupled  differential  equations  can  be  solved  to  determine  the  transmission  of  radia¬ 
tion  using  the  boundary  conditions  /JfO)  =  0  at  the  start  of  the  transmission  path  and  [^{t) 
=  0  at  the  end  of  the  transmission  path. 

A  second  set  of  two-stream  approximations  is  derived  from  eq  9a  and  9b,  which  are  iden¬ 
tical  to  those  presented  by  Bohren  (1987).  This  is  done  by  suppressing  the  term  containing 
the  incident  flux  7;  so  that  the  equations  become 

^/tlr)  =  -/,lr)-^a>.-^i^/r(r)-l-a).-^^^/nr)  (10a) 


^  rw_v  _  r-Yr^  (1  '^^)  r-y  \  (1  5)  r— ^  , 

^  1{\T)  -  /,lr)-wo  — 2 — tiAV-oit — 2 — 


(10b) 


where  the  intensity  /,  now  refers  to  the  total  intensity  (/j  +  1^)  and  not  just  the  diffuse  intensi¬ 
ty.  Equations  10a  and  10b  are  solved  (Appendix  A)  with  the  boundary  conditions  IfiO)  =  I, 
and  /(“(t)  =  0  at  the  start  and  end  of  the  transmission  path,  respectively,  to  yield  the  follow¬ 
ing  expression  for  transmittance  T(T,g) 


T(T,g) 


1 


1  + 


(1-g) 

2 


T 


(11) 


This  equation  is  used  to  describe  the  light  transmission  measurements  that  were  made  in  fall¬ 
ing  snow. 


TRANSMISSION  THROUGH  SNOWFALL 

Visible  light  transmission  data  were  obtained  during  snowfalls  over  a  path  distance  of  650 
m.  1  ne  source  was  a  nearly  parallel  beam  of  light  with  a  diameter  of  0.6  m  and  a  divergence 
of  20  mr.  The  receiver  consisted  of  a  detector  located  at  the  focal  point  of  a  0.6-m-diameter 
lens  with  a  field  of  view  of  0.6  mr.  During  the  transmission  measurements,  the  snow  mass 
concentration  (snow  mass  per  unit  volume  of  air)  was  continuously  measured  at  the  mid¬ 
point  of  the  transmission  path. 

Equation  1 1  shows  that  the  transmission  of  visible  light  through  falling  snow  can  be  deter¬ 
mined  if  the  optical  path  and  the  asymmetry  parameter  are  known.  The  optical  depth  of 
snow-filled  atmosphere  is 

r  =  (12) 

where  A/j  is  the  snow  mass  concentration,  and  the  snow  mass  extinction  coefficient,  a^,  is 
defined  as  the  extinction  cross  section  per  unit  particle  mass.  The  snow  mass  extinction  coef¬ 
ficient  is  typically  around  0.03  mVg  but  it  does  vary  depending  on  the  snow  crystal  type. 

The  asymmetry  parameter  of  snow  was  estimated  to  be  0.9  for  this  study.  The  value  of  0.9 
corresponds  to  the  asymmetry  parameter  in  the  visible  wavelength  for  an  ice  sphere  with  size 
similar  to  a  snow  particle.  Although  the  spherical  approximation  for  snow  is  an  extreme 
idealization,  it  is  used  to  obtain  insights  into  the  scattering  process.  Until  the  phase  function 
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of  snow  is  known,  the  use  of  a  spherical  equivalent  phase  function  may  be  the  most  appro¬ 
priate  approximation. 

From  eq  2,  4  and  1 1 ,  the  transmitted  diffuse  intensity  can  be  expressed  as 


=  - n^r;. - /icxp(-r).  (13) 

1+ 

The  reduced  intensity  (eq  2)  and  the  diffuse  intensity  components  (eq  13)  of  the  transmitted 
intensity  are  calculated  as  a  function  of  snow  mass  concentration  for  a  transmission  path  of 
650  m.  The  results  that  were  calculated  for  a  unit  incident  intensity  are  shown  in  Figure  2.  It 
is  seen  that  the  diffuse  intensity  becomes  dominant  when  the  snow  mass  concentration  ex¬ 
ceeds  0.2  g/mS  which  is  equivalent  to  an  optical  depth  of  3.9. 


Figure  2.  Relative  contributions  of  the  re¬ 
duced  intensity  and  the  diffuse  intensity 
to  the  total  intensity  as  a  function  of  snow 
mass  concentration  for  a  path  distance  of 
650  m. 


Thus  far  it  has  been  assumed  that  all  the  diffuse  radiation  transmitted  in  the  forward  di¬ 
rection  contributes  to  the  total  intensity.  However,  because  of  the  finite  size  of  detectors 
used  to  measure  the  transmitted  light  and  the  spatial  and  angular  spreading  of  the  diffuse 
beam,  only  a  fraction  of  the  diffuse  radiation  is  actually  recorded.  Therefore,  in  order  to 
correctly  interpret  the  experimental  results,  the  measured  diffuse  radiation  /j''  (r.g)  is  ex¬ 
pressed  as 


/r  (’-.«)  =  C(r,g)-/d-(r,g)  (14) 

where  the  correction,  C(T,g),  is  the  fraction  of  diffuse  radiation  that  is  actually  measured. 
The  measured  transmittance,  T‘ir,g),  is  then  expressed  as 


T  (T,g)  =  exp(  -  r)  -t-  C(r,g)  •  /j  (r,g).  (15) 

The  correction  factor  is  dependent  on  the  ratio  of  the  effective  detector  area  to  the  cross- 
sectional  area  of  the  transmitted  diffuse  beam.  In  other  words,  for  a  given  detector  size,  C(T,g) 
decreases  as  the  diffuse  beam  spreading  increases.  Therefore,  information  about  diffuse 
beam  spreading  is  required  to  determine  the  correction  factor.  The  spreading  of  the  diffuse 
beam  is  a  complicated  function  of  the  optical  depth  and  the  snow  phase  function,  therefore 
simplifying  assumptions  are  made  to  determine  the  correction  factor. 

The  diffuse  beam  diameter  increases  as  the  optical  depth  increases.  If  this  increase  is  linear 
with  respect  to  the  optical  dej)th,  C(r,g)  will  be  inversely  proportional  to  the  square  of  the 
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optical  depth.  The  rate  of  increase  is  determined  by  the  phase  function;  the  beam  spreading 
is  less  when  the  phase  function  is  peaked  in  the  forward  direction.  In  other  words,  for  a  given 
optima,  depth,  C(r,g)  will  increase  as  the  asymmetry  parameter  increases. 

.>\nother  consideration  for  estimating  C(T,g)  is  the  single  scattering  limit  at  low  optical 
depths.  At  low  optical  depths,  the  scattered  radiation  that  is  measured  is  dominated  by  a 
sharp  diffraction  peak  in  the  near-forward  direction.  Fraunhofer  diffraction  theory  is  com¬ 
monly  used  to  correct  for  the  effects  of  single  scattering  on  measured  transmission  (Sea- 
graves  1983,  Bohren  and  Koh  1985).  These  calculations  indicate  that  for  visible  radiation,  a 
significant  portion  of  the  light  contained  in  the  first  diffraction  lobe  is  measured  by  a  typical 
transmission  system  (approximately  35  to  40%  of  all  the  scattered  radiation  is  measured  at 
low  optical  depths).  From  eq  8b  it  is  seen  that  a  particle  with  an  asymmetry  parameter  of  0.9 
scatters  95%  of  the  light  in  the  forward  hemisphere.  Therefore,  the  limiting  value  of  C(T,g) 
at  low  optical  depth  (ratio  of  the  scattered  radiation  that  is  measured  to  the  radiation  scat¬ 
tered  in  the  forward  direction)  is  expected  to  range  from  0.37  to  0.42.  A  correction  factor  of 
0.4  is  used  for  the  low  optical  depth  limit  in  this  report. 

By  combination  of  the  above-mentioned  criteria  required  for  the  correction,  the  following 
expression  is  derived  for  the  C(r,g) 


C(T,g) 


0.40 

I  +(1  -g)T^  ■ 


(16) 


Figure  3  illustrates  values  for  C(T,g)  for  several  asymmetry  parameters.  This  figure  shows 
that  eq  16  satisfies  the  criteria  discussed  above.  The  values  for  C(T,g)  decrease  as  the  optical 
depth  increases,  increase  as  the  asymmetry  parameter  increases,  and  also  approach  the  single 
scattering  correction  limit  of  0.40  at  the  low  optical  depths. 

The  utility  of  any  model  is  judged  by  its  ability  to  predict  experimental  results.  Figure  4 
shows  the  comparison  of  the  light  transmission  measured  through  falling  snow  with  those 
calculated  using  the  two-stream  approximation  with  the  correction  factor  applied.  The  figure 
shows  that  the  experimental  results  can  be  described  with  reasonable  acuracy  without  the  use 
of  an  empricial  adjustment.  The  transmission  data  presented  in  Figure  4  were  selected  so  that 
the  snow  crystal  type  and  size  were  similar  through  the  range  of  the  snow  mass  concentration. 


Figure  3.  Correction  for  the  fraction  of 
diffuse  intensity  that  is  measured  by  a 
detector  (eq  16).  The  asymmetry  param¬ 
eters  (g)  are  0.9,  0. 7,  0.5  and  0. 1. 


Figure  4.  Comparison  of  the  transmit¬ 
tance  calculated  from  the  two-stream  ap¬ 
proximation  with  the  experimental  re¬ 
sults.  The  solid  line  represents  the  two- 
stream  results  using  eq  15  and  16. 
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It  is  emphasized  that  the  values  for  C{T,g)  were  determined  using  qualitative  arguments. 
The  good  agreement  with  the  experimental  results  obtained  using  the  qualitative  arguments 
may  be  coincidental.  More  work  is  required  to  justify  the  analytical  expression  for  Cfr.g). 
Nevertheless,  it  is  important  to  point  out  that  a  correction  factor  that  is  qualitatively  similar 
to  the  one  presented  is  required  to  interpret  light  transmission  measurements  through  an  op¬ 
tically  thick  snowfall. 


CONCLUSIONS 

The  parameters  (the  optical  depth  and  the  phase  function)  that  are  required  for  calculating 
the  radiative  properties  of  falling  snow  can  only  be  approximated.  Therefore,  complex  and 
time-consuming  methods  to  solve  the  exact  radiative  transfer  equations  are  not  necessary 
when  one  is  interested  in  interpreting  experimental  results.  A  simple  two-stream  approxima¬ 
tion  to  radiative  transfer  where  the  light  is  constrained  to  travel  in  two  directions  only  (for¬ 
ward  and  backward)  is  used  to  describe  the  transmission  of  visible  light.  A  correction  factor 
that  accounts  for  the  portion  of  the  diffuse  radiation  that  is  actually  measured  by  a  detector 
is  then  applied.  The  comparison  of  the  light  transmission  data  obtained  in  snowfall  with 
those  calculated  using  the  two-stream  approximations  is  in  good  agreement  over  a  wide 
range  of  optical  depths. 
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APPENDIX  A:  TWO-STREAM  APPROXIMATION  SOLUTION 


Adding  and  subtracting  the  following  differential  equations 

d  N  ,  (1  +«)  ,  ,  X  ,  (1  -«)  X 

^/(  (t)  =  -/[  (t)  +  u'o  — 2 — ■'(  +  — 2 — 

d  .  (1  +g)  ,  (1  -g)  . 

^  /[(’■)  =  V,  (t)-cu'o  — 2 — A  — 2 — 

result  in  the  following; 

(/f+Zf)  =  (1 -a.-o^)(/r-/r) 

^  (zr-zr)  =  (-1 +u.’o)(zr+zr). 

For  the  case  where  absorption  is  neglible  (a-o  =  1),  eq  A2a  and  A2b  become 

^  (Zf+  zf )  =  (1  -g)ur-  in 

i  A-)  -  0. 

The  general  solutions  to  eq  A3a  and  A3b  are 

zr-zr=  c, 

Zt-+Zf=  G-(l  -g)-CrT 
which  can  be  reexpressed  as 

Zf=  y(C,+C)-  -C.-r 

Zr=  y(C,-Ci)-  C.  r 

where  constants  C,  and  C:  are  determined  by  the  boundary  conditions. 
Using  the  boundary  conditions 

1\(t)  =  Zj  for  t  =  0 

Zt'tr)  =  0  for  T  =  (7„,  •  L 
the  following  is  derived: 

C|  +  C2  =  2  Zj 

C:-C,  =  (l-g)-C,-r. 
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(Ala) 

(Alb) 

(A2a) 

(A2b) 

(A3a) 

(A3b) 

(A4a) 

(A4b) 

(A5a) 

(A6b) 


(A7a) 

(A7b) 


These  equations  can  be  manipulated  to  show  that 

h 


c,  = 


(A8) 


The  transmitted  radiation,  T,  is  then  expressed  as 
/f  1 


T  = 


(A9) 
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